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ABSTRACT 


In  this  thesis,  a  method  for  the  evaluation  of  the  transient 
response  of  a  nonuniform  lossy  transmission  line  is  presented.  The 
method  of  characteristics  is  employed  to  obtain  difference  equations 
describing  the  transmission  line. 

A  stepped  line  approximation  is  used  to  analyze  the  tran¬ 
sient  response  of  the  given  nonuniform  line.  The  concept  of  electrical 
length  is  employed  in  dividing  the  line  into  a  number  of  equal  delay 
sections.  The  set  of  difference  equations  describing  the  set  stepped 
line  is  suitable  for  digital  computer  solution. 

The  computational  procedure  involving  the  use  of  a  digital 
computer  is  illustrated  for  a  specific  distributions  of  L(x),  C(x), 
r (x)  and  g(x). 
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CHAPTER-I 


Introduction 

1.1  Background 

In  the  design  of  all  apparatus  and  machines  related  to 
the  generation  and  transmission  of  electrical  power,  it  is 
important  to  know  the  salient  features  of  the  overvoltages  likely 
to  be  encountered.  Faults  caused  by  switching  are  becoming  more 
important  with  the  ever  higher  transmission  voltages  being  used. 

Further  the  basic  form  of  implementing  computer  hard¬ 
ware  is  in  mounting  integrated  circuits  (I.C.)  packages  on 
multilayer  printed  circuit  boards  which  in  turn  are  mounted  on 
so  called  mother  boards.  The  pins  of  the  I.C.  packages  on  the 
printed  circuit  board  are  usually  connected  to  each  other  by  etched 
metal  strips,  and  circuit  engineers  are  forced  to  regard  these 
metal  strips  as  a  length  of  a  transmission  line  as  the  speed  of 

(g) 

the  circuit  approaches  the  nanoseconds  range v 

Furthermore,  the  problem  of  constructing  a  mathematical 
model  for  simulating  the  electrical  properties  of  memory  arrays 
has  been  considered  in  1963  by  Weeks,  who  found  that  a  memory 

array  can  be  approximated  by  a  system  of  linear  transmission  lines. 

Because  of  the  above  considerations,  transient  analysis 
of  general  transmission  lines  has  recently  become  the  object  of 


1 


2 


wide-spread  interest  of  many  investigators.  The  analysis  of  trans¬ 
ient  waves  on  a  transmission  line  leads  to  a  set  of  partial 
differential  equations  which  define  the  relation  between  current 
and  voltage  as  functions  of  time  and  position  on  the  line.  In 
general,  these  equations  do  not  lend  themselves  to  a  closed  form 
solution  except  in  the  special  case  of  distortionless  lines  (where 


£ 

—  =  G/C),  which  includes  the  lossless  case  (R  =  G  =  0).  Thus  using 
numerical  techniques  for  solving  such  equations  is  inevitable 
The  application  of  the  method  of  characteristics  was 


shown  in  1967  by  F.H.  Branin  Jr., 


(1) 


to  provide  a  simple  analytic 


solution  for  the  problem  of  a  uniform  lossless  line.  The  method 

was  shown  to  be  superior  in  both  speed  and  accuracy  to  the  more 

familiar  method  of  integrating  the  differential  equations  that 

describe  a  lumped  L.C.  model  of  the  line.  Branin  did  not  attempt 

in  this  work  to  extend  the  method  to  more  general  lines. 

(2) 

Y.K.  Liu,  in  1968,  considered  the  transient  analysis 

of  a  uniform  lossy  line  terminated  in  a  resistive  load.  For  his 

analysis  he  used  the  method  of  characteristics  and  a  second  order 

Runge-Kutta  technique  to  solve  numerically  the  resulting  pair  of 

ordinary  differential  equations. 

Transform  methods  for  analyzing  the  transient  response 

(3) 

of  uniform  lines  was  considered  by  K.A.  Chen  in  the  same  year. 

Chen's  work  was  concerned  with  interconnected  lines  as 


an  aid  in  the  design  of  memory  arrays  in  digital  computer  hardware. 


3 


(4) 

The  investigation  carried  out  by  Wassel  used  transform 
methods  for  the  analysis  which  was  concerned  with  the  effects  of  RC 
loadings  of  pulse  signal  transmission  lines.  It  is  noted  that 
R.  Murray- Shelley investigated  the  same  problem  in  the  same 
year  using  the  method  of  characteristics  (often  called  the  graph¬ 
ical  method).  In  1969,  Y.K.  Liu^^  reported  the  application  of 
the  same  technique  he  used  in  (2)  when  the  line  is  terminated  in 


a  tunnel  diode. 


(9) 


In  1970,  V.  Dvorak  reported  a  novel  method  of  treating 
this  problem.  He  used  the  method  of  characteristics  in  a  way  that 
is  straight  forward  and  avoids  entirely  lengthy  numerical  techniques 
such  as  Runge-Kutta  method.  One  of  the  advantages  of  Dvorak's 
method  is  that  any  loading  configuration  can  be  handled  easily. 
Further  the  treatment  of  nonuniform  lines  was  shown  to  be  best 
analyzed  using  a  stepped  line  approximating  the  original  line, 
that  is  the  division  of  the  line  into  a  number  of  sections  having 
the  same  delay  (electrical  length) . 

1.2  Scope  of  the  thesis 


In  this  thesis  the  application  of  the  method  of  characteristics 
to  the  problem  of  computing  the  transient  response  of  a  general 
transmission  line  is  presented.  The  analysis  for  a  nonuniform  loss¬ 
less  line  is  presented.  The  method  of  analyzing  both  uniform  and 
nonuniform  lossy  lines  is  discussed,  here  a  stepped  line  approx¬ 
imating  the  given  line  is  used.  The  choice  of  proper  time  step  and 


4 


number  of  sections  of  the  stepped  line  is  considered  for  the  case  of 

uniform  lines.  This  is  to  some  extent  an  extension  of  the  work  of 

(9) 

V.  Dvorak. 

A  general  transmission  line  transient  analysis  computer 
program  is  described.  The  program  is  written  in  FORTRAN  IV  pro¬ 
gramming  language.  The  program  is  discussed  in  some  detail  and 
numerical  results  are  given  for  some  example  transmission  lines. 


CHAPTER  2 


Transient  analysis  of  lossless 


nonuniform  transmission  lines 


2.1  The  method  of  characteristics 

The  set  of  partial  differential  equations  describing  a 
uniform  transmission  line  is, 


3V 

3Z 


(2.1) 


(2.2) 


C  9V  +  G  V 

at 


where  L,  C,  R  and  G  are  inductance,  capacitance,  resistance  and 
conductance  per  unit  length  of  the  line  respectively,  V  =  V(Z,t) 
and  i  =  i(Z,t)  are  the  voltage  and  current,  respectively,  at 
distance  Z  from  one  end  of  the  transmission  line  at  time  t.  In 
essence  the  method  of  characteristics  transforms  (2.1)  and  (2.2) 
into  two  ordinary  differential  equations  as  follows.  Let 


Ri  +  Li  +  V 


Z 


0 


(2.3) 


GV  +  CV  +  i 

L  Lj 


0 


(2.4) 
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where 


3V 

3Z 


V 


Z 


V 

t 


3V 
3 1 


3_i 

3Z 


i 


t 


3  i 
3 1 


Def ining 

4>  =  <j>1  +  X  4>2  =  0 


Then 


A  V7 

<f>  =  Ri  +  A  GV  +  L[it+^iz]  +  A  C[Vt-t~] 


=  0 


(2.5) 


This  is  a  voltage  equation  so  that  the  quantity  between 
brackets  in  the  third  term  of  (2.5)  should  be  equal  to  the  total 
derivative  of  the  current  3^,  hence  we  have 


di 

dt 


but 


di  =  .  3i  dZ 

dt  Xt  3Z  dt 


hence 


7 


Similarly  the  quantity  between  brackets  in  the  last 

term  of  (2.5)  should  be  equal  to  the  total  derivative  of  the 
dV 

voltage  —  ,  hence 
dt 

d_Z  _  1_ 

dt  XC  (2.7) 


by  (2.6)  and  (2.7)  we  have 


(2.8) 


or 


dZ 

dt 


Substituting  (2.8)  into  (2.5) 


di  .  rr ~  dV 


4>  =  Ri  +  GV  y  —  +L-1VLC 


=  0 


or 


k  'v  Vt 


“  [77  V  ±  —  i] 

c  /lc 


(2.9) 


(2.10) 


Now  (2.9)  implies  that  in  the  (Z-t)  plane  we  have  essentially 
straight  line  characteristic  curves  in  the  case  of  a  uniform  trans¬ 
mission  line. 


202  The  stepped  line  approximation 

Consider  a  lossless  transmission  line  characterized  by 
distributed  inductance  L(Z)  >  0  and  capacitance  C(Z)  >  0  per  unit 
length,  where  Z  is  the  physical  position  on  the  line.  The  electrical 


' 


8 


position  along  the  line  is  defined  by: 

Z 


y(Z)  = 


I  < 
0 


l  (n)  c  (n)  d, 


(2.11) 


and  the  local  characteristic  impedance  by 
P  (y)  =  VL  [Z  (y)  ]  /  C  [Z  (y)  ] 


(2.12) 


Now  the  characteristic  curves  in  the  (Z-t)  plane  are  not 
straight  lines  in  the  case  of  nonuniform  transmission  lines,  but 
the  introduction  of  the  electrical  position  y(Z)  yields  straight 
line  characteristic  curves  in  the  (y-t)  plane. 


A  transmission  line  of  total  length  h  may  be  divided 

If  n  is  chosen 

large  enough,  then  any  such  section  can  be  approximated  by  a 


into  n  sections  of  the  same  delay  Ay  =  ^ 

n 


uniform  lossless  line  with  an  average  characteristic  impedance 
p  (1  ^  K  $  n)  and  with  total  propagation  delay  T  =  Ay .  The  result- 
ing  cascade  of  n  uniform  lossless  terminating  networks  is  described 
by  2  n  difference  equations  in  the  time  domain  of  the  form  given 
in  (A-9)  and  (A-10),  plus  the  two  equations  describing  the  termin¬ 
ating  networks.  A  stepped  line  approximating  the  given  nonuniform 
line  is  shown  in  figure  2.1. 

The  basic  equations  for  a  uniform  lossless  line  are  derived 
in  Appendix  A.  Applying  (A-9)  to  the  incremental  line  between  the 
(K-l)St  and  nodes  in  figure  2.1  yields 

MO  +  PK_X  iK(t)  =  VK_1  (t-T)  +  pK_1  iK_1  (t-T)  (2.13) 


K-l 


2^K$n+l 


I  V  I 


9 


where  we  substituted 


VR(t)  =  V  (1 ,  t) 


VK-1  (t)  =  V(0’t) 


Further  (A-10)  is  applied  to  the  incremental  line  be¬ 
tween  the  Kt^1  and  (K-l)St  nodes  giving 

VK('t')  “  PK  1K('t')  =  VK+1  ('t_T^>  ”  PK  XK+1  ('t_T')  (2.14) 

l.$K^n 


The  terminating  networks  A  and  B  are  given  by  their  V-A 
characteristics 


A  :  i_ 


f  (V  ) 

i  v  r 


B  :  i  , .  =  f  (V  ) 

n+1  2  n+1 


(2.15) 

(2.16) 


Given  the  initial  conditions  at  the  time  instant  t=0, 

the  values  of  V  (t)  and  i  (t)  can  then  be  computed  at  time 

K  K 


instants  t  =  mT,  l<:m$m 

max 
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PIG .  2.1  A  STEPPED  LINE 


CHAPTER  3 


Transient  analysis  of  lossy  uniform  lines 
3.1  Analysis 

Consider  a  lossy  uniform  line  for  which  (2.9)  implies 
that  in  the  (Z-t)  plane  we  have  essentially  straight  line 
characteristic  curves  as  shown  in  figure  (3.1). 


where  along  the  a  curve  we  have 
dZ  1 

dt  =  VIZ 

while  along  the  3  curve  we  have 


dZ^  =  1_ 

dt 

Using  (2.10)  we  have  the  following  two  differential  equations  each 
being  defined  along  a  different  characteristic  direction  in  the 
(Z-t)  plane 

A1  dZ  1 

Along  the  a  curve:  - —  =  — 

dt  /lc 

4t  [v+Pi]  =  -  1=  [pGV+Ri]  (3.1) 

dt  JTc 


Along  the  6  curve  :  d_Z 

dt 


d_ 

dt 


[V-pi] 


[pGV-Ri ] 


(3.2) 
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Note  that  (3.1)  defines  a  forward  propagating  wave  where  as  (3.2) 
defines  a  backward  propagating  wave. 

Dividing  the  line  into  n  sections  of  delay 
T  =  h  v|lcT  /  n  (3.3) 

and  taking  dt  =  At  =  T,  one  gets 

K+i 

A[V  (Z,t)  +  pi  (Z ,  t )  ]  =  —  [G  p  V  (Z,t) 

n 

K 


+  R  i  (Z,t)] 

K  -h 

A [V  (Z,t)  -  pi  (Z,t)]  =  —  [G  p  V  (Z,t) 

K+l  n 

-  R  i  (Z,t)] 


(3.4) 


(3.5) 


FIG.  3.1  FAMILIES  OF  a  AND  g  CHARACTERISTIC 


CURVES  IN  Z-t  PLANE 


letting 
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Gh  Rh 

g  =  —  ,  r  =  — 
n  n 


(3.6) 


Then  (3.4)  and  (3.5)  become 


K+l 

A  [V  (Z  ,  t )  +  p  i  (Z ,  t )  ]  =  -  [g  P  V  (Z,t) 

K 

+  r  i  (Z  ,  t )  ] 


(3.7) 


K 


A  [V  (Z ,  t )  -  p  i  (Z,t)l  =  -  [g  p  V  (Z,t) 

K+l 


-  r  i  (Z ,  t )  ] 


th 


(3.8) 


Notice  that  any  event  occuring  at  the  K  node  at  the 

s  t 

time  instant  (t-T)  arrives  at  the  (K+l)  node  at  the  time  instant, 
t,  for  the  forward  propagating  wave,  (the  reverse  is  true  for  back¬ 
ward  propagating  wave). 

So  that  writing  (3.7)  and  (3.8)  in  a  discrete  (difference 
equation)  form  leads  to  the  following  question: 

how  are  we  going  to  discretize  the  right  hand  side  of 
each  of  (3.7)  and  (3.8)?  To  this  end  we  let 


f  (Z ,  t )  =  g  p  V  (Z  ,  t )  +  r  i  (Z,t) 


(3.9) 


g (Z  ,  t )  =  g  p  V  (Z,t)  -  r  i (Z , t ) 


(3.10) 


We  have  the  following  two  cases 

I.  If  we  assume  that  f(Z,t)  and  g(Z,t)  remains  constant  over  the 
time  interval  [t-T,  t],  this  is  equivalent  to  the  assumption 


f  (Z  ,  t )  =  fK(t-T) 


(3.11) 
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g(z>t)  =  gK+1(t-T) 


(3.12) 


II.  If  we  assume  that  f(Z,t)  and  g(Z,t)  change  in  a  linear  fashion 


over  the  time  interval 

f 

f(Z,t)  =  fR(t-T)  +  - 


[t-T,t]  then 

K+l  ( t )  ~  fK(‘t~T') 

2 


or 


f(Z,t)  =  [ f K (t-T)  +  fK+1(t)]  /  2 


(3.13) 


and 


g(Z,t)  =  [gR+1(t-T)  +  gK(t)]  /  2  (3.14) 

Application  of  (3.11)  and  (3.12)  to  (3.7)  and  (3.8)  respectively 
yields 

VR+1(t)  +  p  iK+1(t)  =  [1  -  Pg]  VK(t-T)  + 

[p  -  r]  iK(t-T)  (3.15) 


VR(t)  -  p  iR(t)  =  [1  -  Pg]  vK+1(t-T)  - 

[p  -  r]  iK+1(t-T)  (3.16) 

Equations  (3.15)  and  (3.16)will  be  referred  to  as  the  first  form¬ 
ulation  equations. 

Further,  application  of  (3.13)  and  (3.14)  to  (3.7)  and 


(3.8)  respectively  yields 


K  * 


15 

VK+I(t>  '1+£f]  +  Wc)  [p  +f]  ■ 

VK(t-T)  [1  -  £|]  +  iK(t-T)  [p  -  |]  (3.17) 


VR(t)  [1  +  £§]  -  lR(t)  [p  +  f]  = 

VK+l(t'T)  l* 1  -  -  iK+l(t'T)  [P  *  1]  (3.18) 

Equations  (3.17)  and  (3.18)will  be  referred  to  as  the  second  form¬ 
ulation  equations. 

The  first  and  second  formulation  equations  are  the 
difference  equations  of  a  very  short  lossy  line  in  the  time  domain. 

3.2  Boundary  conditions 

Let  the  line  have  the  following  initial  conditions 
iR(0)  =  0  UK<n  +  1  (3.19) 

VT,(0)  =  0  l^K^n  +  1  (3.20) 

K 


Let  the  sending  end  terminating  network  be  an  ideal  current  source 

whose  current  output  is  a  prespecified  function  of  time  I(t).  If 

the  source  is  assumed  to  have  a  shunt  conductance  G  ,  then  the 

s 

following  equation  holds  true. 


i  (t)  =  I(t)  -  G  V  (t) 

1  si 


(3.21) 


If  the  line  is  terminated  in  a  load  resistance  R  then 

r 


Vnii)  -  Rr  W 


(3.22) 


The  line  with  its  terminating  network  is  shown  in  figure  (3.2). 
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FIG.  3.2  UNIFORM  TRANSMISSION  LINE  WITH 

ITS  TERMINATING  NETWORKS 
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3.3  Implementing  the  first  formulation 

We  rewrite  the  first  formulation  equations  (3.15)and 

(3.16)  as 

VK(t)  +  P  iK(t)  =  [1  -  pg]  V^ft-T]  + 

[p  -  r]  iK_1(t-T)  (3.23) 


vK(t)  “  P  iK(t)  =  [1  -  Pg]  vK+1(t-T)  - 

[p  -  r]  iK+1(t-T)  (3.24) 

Adding  (3.24)  to  (3.23)  and  subtracting  (3.24)  from  (3.23)  re¬ 
spectively  yields 

vK(t)  =  (1  -~--P£)  [VK_1(t-T)  +VK+1(t-T)]  + 

[iR_1(t-T)  -  iK+1(t-T)]  2<R$  (n-1)  (3.25) 

-  (Hr&>  [vk-i(c't)  -  vk+i(c't)1  + 

h[l  ~  f]  [iK_1(t-T)  +  iK+1(t-T)]  2<K$  (n-1)  (3.26) 

Thus  the  values  of  VT.  (t)  and  iTr(t)  at  times  t  =  mT ,  l<:m$m  can 

K  K  max 

be  computed  using  (3.23),  (3.24),  (3. 25)4(3.26) . 


■ 
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3.4  Implementing  the  second  formulation 


Rewriting  the  second  formulation  equations  (3.17)  and 

(3.18)  as 

VR(t)  [1  +  ^f]  +  iK(t)  [p  +  §]  =  VK_1(t-T)[l  -  £|] 


+  iR_1(t-T)  [p  -  |] 


(3.27) 


VK(t)  ^  +^]  ~  iK(t)  tp  +  |l  =  VK+i(t"T)[1  -  £f] 


1K+l('t  T')  2  ^ 


(3.28) 


Adding  (3.27)  and  (3.28)  and  subtracting  (3.28)  from 
(3.27)  respectively  yields 


vK(t)  =([vK_1(t-T)  +  VK+1(t-T)]  [1  -  ^f]  + 
[iK_1(t-T)  -  iR+1(t-T)]  [p  -  | ]^  /  (2  +  pg) 

={[vK-i(t-T)  -Vi(t-T)]  [1  + 

[iK_1(t-T)  +  iK+1(t-T)]  [p  -  j  (2p  +  r) 


(3.29) 


(3.30) 


Thus  the  values  of  VR(t)  and  iK(t)  at  instants  t  =  mT ,  Um^m^^ 
can  be  computed  using  equations  (3.19),  (3.20),  (3.29),  (3.30), 


(3.21),  (3.22). 


3.5  Error  estimation 


If  we  let 

W(Z,t)  =  V(Z,t)  +  p  i(Z,t) 


(3.31) 


f(Z,t)  =  —  [G  p  V(Z,t)  +  R  i(Z,t)]  (3.32) 

J  LC 

Then  equation  (3.1)  is  transformed  to 


dW 

dt 


=  f(t) 


(3.33) 


for  a  specific  position  Z. 

Thus  the  increment  of  change  of  W,  which  is  denoted  by 

AW  over  the  interval  [t  ,  t  +  At ]  is 

o  o 


AW 


[f  (t)  ] 


avg 


Where  [f  (t) ] 

interval  [t  ,  t  +  At] 
o  o 


[f (t)  ] 


avg 


1_ 

At 


.  At  (3.34) 

is  the  average  value  of  f(t)  in  the 

avg 

given  by 

t  +  At 
C  ° 

f(§)  d  §  (3.35) 


t 

o 


Now  using  Taylor’s  expansion  for  the  right  hand  side  of 


(3.35)  one  gets 
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[f (t)  ] 


avg 


£(to> 


(At)2 

3! 


f(t)  |t  + 
o 


f(t) 


.  + 


(At) 


.n 


(n+1) 


(3.36) 


If  At  is  sufficiently  small  then 


[f (t)  ] 


avg 


f(to)  + 


At 

2! 


f(t) 


o 


+ 


(At)2  " 
3! 


(t) 


t 

o 


+ 


(At)3 

4! 


T(t) 


t 

o 


(3.37) 


Thus  a  sufficiently  accurate  value  of  AW  is 


AW  =  A [ f  +  4 
o  2 


f 

o 


0  f* 


(3.38) 


where 


f 

o 


f 

o 


f  (t) 


t 

o 


and  so  on. 


a  -  First  formulation 


Referring  to  (3.11),  the  average  value  of  f(t)  over 

[t  ,  t  +A]  was  taken  to  be  f(t  ),  (where  t  =t-T,  A=T  in  this 
o  o  o  o 

thus  the  estimated  increment  of  change  in  W,  denoted  by  AW^ 


the  interval 

case) 

is 


[f  (t)  ] 


avgx 


.  At 


(3.39) 
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hence 

AW  =  f  .  A  (3.40) 

o 

This  provides  us  with  a  sufficiently  accurate  error 
measure  for  the  first  formulation  which  we  denote  by  e  ,  so  we 

S  1 

have 


e  =  AW  -  AW  (3.41) 

S1  S1 

Substituting  (3.38)  and  (3.40)  into  (3.41)  we  get 
A2  ‘ 

e  =  ~  f  +  0,  (At)  (3.42) 

S.,  Z  O  J 


with  O^(At)  being  terms  in  At  higher  than  third  order. 

Now  if  an  upper  limit  on  the  error  incurred  in  im¬ 
plementing  the  first  formulation  is  given  by  e  ,  then, 

a 


£ 

a 


f 

o 


or 


At  < 


f 

o 


(3.43) 


which  defines  an  upper  limit  on  the  value  of  the  time  step  to  be 
chosen  so  that  the  given  accuracy  is  achieved.  Let  this  upper  limit 
be  At-^,  then 


o 
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b  -  Second  formulation 


Referring  to  (3.13),  the  average  value  of  f(t)  over  the 
interval  [tQ,  t  +A]  was  taken  to  be 


[f(C)1avg2  =  l2[£o  +  fl] 


(3.44) 


where 


f1  =  f(t  +  At) 
1  o 


(3.45) 


thus  the  estimated  increment  of  change  in  W  denoted  by  AW  is 

S2 


AW  =  [f (t)  ] 


avg. 


At 


(3.46) 


which  yields  the  error  measure  for  the  2nd  formulation  e  which  is 

S2 

given  by 


=  AW  -  AW 
S2  S2 


(3.47) 


applying  (3.38)  and  (3.46)  to  (3.47)  yields 


e  =  ~  f  +  0  (At) 

s^  12  o  4 


(3.48) 


Now  if  an  upper  limit  on  the  error  incurred  in  implementing  the  2nd 
formulation  is  given  by  then 


A  " 

£  >  -  f 

a  12  o 


or 


At 


12e. 


1/3 


(3.49) 


23 


def ine 


At2  = 


12E; 


1/3 


o 


This  specifies  the  upper  limit  (At^)  on  the  value  of  the  time  step 
to  be  chosen,  so  that  the  given  accuracy  is  achieved, 
c  -  Simplified  calculation  of  f 

o 


In  appendix  (B)  it  is  shown  that  the  function  W  obeys  the 


following  equation 


WK+1  =  WR  exp[-At.s-M  —  +  gp )  ] 


(3.50) 


where  s  is  the  Laplace  operator  and  hence  the  term  At . s  in  the  exponent 
represents  the  time  delay  between  two  consecutive  nodes. 

We  are  concerned  only  with  the  deterioration  in  the  wave  W, 


not  the  delay,  so  that  we  need  only  to  consider 


W, 


K+l 


W, 


J  exp [-^(  ^  +  gp)] 
R  P 


(3.51) 


Let 


e  =  exp  [-%  (  ■—  +  gp)] 


(3.52) 


Then 


W  =  e  W 
K+l  K 


(3.53) 


2 

W  =  e  W 
K+2  K 


(3.54) 


WK+3  ‘  ^  WK 


(3.55) 


24 


This  enables  one  to  evaluate  f  and  f  which  are  required  in  the 

o  o 

calculation  of  upper  limits  on  the  time  step  in  the  previous  section. 


By  definition 


r/  N  dW 
f(t)  =  dF 


(3.56) 


Then  in  a  discrete  form  one  can  write 


f  = 
o 


W  -  W 
K+l  K 

At 


(3.57) 


Substituting  (3.53)  into  (3.57)  then 


W 

K 

fc=AF  t*'1’ 


(3.58) 


where 


f  =  f(o) 
o 


(3.59) 


and 


h  - 


W  -  W 
K+2  K+l 

At 


(3.60) 


Substituting  (3.53)  and  (3.54)  into  (3.60)  one  gets 


WK 

£i  -  ZF  £(E-X) 


(3.61) 


Similarly 

W. 


f  2 


K 

At 


£  (e-1) 


where 


f  =  f(m.At) 
m 


(3.62) 


(3.63) 


further 
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f(t) 


d  f  (t) 

dt 


in  a  discrete  form 


fr£ 


o 


At 


hence  using  (3.58)  and  (3.61) 


f 

o 


W 


K 


(At) 


(e-1) 


2 


and 


f  Wr  <  ^^2 

f  =  - x  e(e-l) 

1  (At)Z 


where 


f 

m 


f  (t) 


(m. At) 


since 


f  (t) 


d  f  (t) 
dt 


(3.64) 


(3.65) 


(3.66) 


(3.67) 


(3.68) 


(3.69) 


in  a  discrete  form 


f 

o 


£r£ 


o 


At 


hence,  using  (3.66)  and  (3.67)  one  gets 


f 

o 


W 


K 


(At) 


(e-1) 


3 


(3.70) 


(3.71) 
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d  -  proper  number  of  sections 

The  results  obtained  in  the  last  two  sections  will  now  be 
applied  to  obtain  an  estimate  of  the  minimum  number  of  sections  to 
which  the  line  should  be  divided  in  order  to  meet  a  specified  accuracy 
limit . 

Let  the  total  allowable  error  be  AW^  defined  as  : 

A  W  =  n  .  e  (3.72) 

t  a 

Denote  the  relative  error  in  calculating  W  by  p  so  that 
AW 

u  =  “j r  (3.73) 

Also  we  have  by  (3.52) 

c  =  exp [-  —  ]  (3.74) 

n 

where 

a  =  -  +  G  p]  (3.75) 

P 

1  ex 

In  the  practical  case  where  r  <<  p  and  g  <<  — ,  —  <<  1 

and  e  can  be  approximated  by 

c  =  1“-  (3.76) 

n 

Using  this  expression  for  e  we  obtain  the  following  results. 

A  -  First  formulation 

Substituting  (3.66)  into  (3.43)  one  obtains 

(1-0  <  (2e  /W  )** 

a  n 


(3.77) 


27 


further  substitute  (3.72)  and  (3.73)  into  (3.77)  to  obtain 


(1-e)  <' 


2jj 


(3.78) 


Using  (3.76)  in  (3.78)  one  obtains 


a 

—  < 

n  Y  n 


2y 


(3.79) 


or 


n  > 


a 

2y 


(3.80) 


Thus  the  minimum  number  of  sections  to  be  chosen  for  the  first  form¬ 
ulation  is  given  by 


a 

nl  "  2y 


(3.81) 
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B  -  Second  formulation 


Substituting  (3.71)  into  (3.49)  one  obtains 


/ 


1-e  < 


12e 


1/3 


W 


(3.82) 


K 


We  obtain  upon  substitution  of  (3.72)  and  (3.73)  into  (3.82)  we  get 

1/3 

1-e  <  —  (3.83) 


n 


Then  (3.76)  and  (3.83)  yield 

1/3 


SL  < 

n 


12p 

n 


(3.84) 


or 


n  >  [a  /12p] 


(3.85) 


Thus  the  minimum  number  of  sections  to  be  chosen  for  the  second 


formulation  is  given  by 


n^  = 


a 


12m 


(3.86) 


3.6  Comparison  between  the  two  formulations 

In  this  section  it  is  shown  that  for  most  practical  trans¬ 
mission  lines  to  achieve  the  same  accuracy  by  using  the  two  given 
formulations,  the  first  formulation  requires  a  larger  number  of 
sections  than  that  required  by  the  second  formulation. 

Let 


(3.87) 
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where  n^  and  are  the  minimum  number  of  sections  for  a  given 
accuracy  measure  y  required  for  the  first  and  second  formulations 
respectively  as  given  by  (3.81)  and  (3.86). 

Hence 


a  r  a  1  . 
p  =  47  [  7  '  3  1 


(3.88) 


The  required  relative  error  normally  does  not  exceed  0.01,  thus  we  can 
assume 


y  $  0.01  (3.89) 

Now  by  definition 

a  =  hi  ~  +  Gp]  (3.90) 

P 

so  that  using  (3.89)  and  (3.90) 

If 

7+Gp>l50  (3-91) 

then 

p  >  0  and  n^  >  n^ 

table  3.1  gives  a  comparison  between  n^  and  n?  for  different  values 


of  a . 
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a 

0. 99999964E-01 
0 . 19999993E  00 
0  29999989E  00 
0.39999986E  00 
0.49999987E  00 
0  59999979E  00 
0.69999975E  00 
0  -  79999971 E  00 
0.89999968F.  00 
0.99999964E  00 


0.49999967E-01 
0.19999987E  00 
0.44999969E  00 
0.79999948E  00 
0.12499990E  01 
0.17999983E  01 
0.24499979E  01 
0.31999979E  01 
0-40499983E  01 
0 . 49999981E  Qi 


0. 28867  505E-01 
0.81649601E-01 
0 . 1 4999998E  00 
0.23094004E  00 
0.32274854E  00 
0.42426401E  00 
0.53463376*  00 
0.65319717E  00 
0.77942276*  00 
0 . 91287082E  00 


TABLE  3.1  VARIATION  OF  n  AND  n2  WITH  a  FOR  y=0.1 


\ 
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CHAPTER  4 


TRANSIENT  ANALYSIS  OF  LOSSY  NONUNIFORM  LINES 


4.1  Dividing  the  line 


Consider  a  lossy  transmission  line  characterized  by  dis¬ 


tributed  inductance  L(Z)  >  0,  capacitance  C(Z)  >  0,  resistance 

R(Z)  >  0  and  conductance  G(Z)  >  0  per  unit  length,  where  Z  is  the 

physical  position  on  the  line.  In  this  case  the  relations  given  by  equations 

(3.1)  and  (3. 2)  hold  true,  but  with  the  transmission  line  parameters 

varying  with  Z.  This  means  that  the  characteristic  curves  in  the 

Z-t  plane  are  no  longer  straight  lines  as  is  the  case  for  uniform  lines. 


The  electrical  position  along  the  line  is  defined  by 


Z 


(4.1) 


0 


From  this  definition  it  is  evident  that  the  characteristic  curves 
in  the  Y-t  plane  are  straight  lines.  It  is  worth  noting  here  that 
the  electrical  position  is  the  time  delay  that  a  wave  initiated  at 
Z=0  takes  to  arrive  at  the  physical  position  Z. 

Thus  if  we  have  a  line  whose  length  is  h,  the  total  time 
delay  would  be 


0 


(4.2) 


Now  if  the  physical  line  length  was  divided  into  any  number 


of  equal  sections,  the  electrical  length  (delay)  of  these  sections  will 
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not  be  equal.  Hence  we  need  to  divide  the  line  into  n  sections  of 
equal  electrical  length,  the  method  that  will  be  used  is  as  follows: 

I  -  Divide  the  physical  line  length  h  into  N  sections, 
so  that  each  section  is  of  physical  length  —  . 

Thus  we  have  created  an  N-component  vector  [Z  ]  representing  the 

K 

s  t 

physical  length  between  the  sending  end  and  the  (K  +  1)  node. 

With  Z  =  (K  -  1) .  —  K  =  2,...,  N+l  (4.3) 

K.  In 

II  -  Using  (4.1),  another  N-component  vector  [Y  ]  representing 
the  electrical  length  between  the  sending  end  and  the  (K+l)st  node 
corresponding  to  the  above  divisions. 

III  -  Now  the  electrical  line  length  £  is  divided  into  n 
equal  sections.  The  delay  of  each  is  At. 


At  =  — 
m 


(4.4) 


Thus  a  new  n-component  vector  [Y^.]  is  created 


At 


j  !»••••»  ^ 


(4.5) 


IV  -  The  physical  position  vector  [Z^ ]  corresponding  to 
the  vector  [Y^ ]  is  obtained  using  the  linear  interpolation  formula 


Vl  -  ZK  + 


<ZK+1  -  V  (Yi  -  V 

(yk+i  -  V 


j  !»•••» 


n-1  (4.6) 


Thus  the  line  is  now  divided  into  n  sections  of  equal 


delay. 

4.2  The  main  equations 

Rewrite  equation  (3.1) 

ft  i(z>t>i  = 


rG(Z) 

LC(Z) 


V(Z,t)  ± 


R(Z) 

v/l( z)  C(Z) 


i  (Z,t)] 


(4.8) 


Let 

W(Z,t)  =  v  (Z  ,t ) 


(4.9) 


f(Z,t)  =  -  [  p&r  V  (Z,t)  +  R(Z)-— .  -  i  (Z ,  t )  ]  (4.10) 

sj L(Z)  C(Z) 


Thus  one  of  the  two  equations  given  by  (4.8)  is  trans¬ 


formed  to 


d  W(Z,t) 

dt 


f(Z,t) 


(4.11) 


Consider  the  time  interval  [t  -  T,  t]  with  T  =  At  being 
the  time  delay  of  each  section  as  given  by  (4.4),  then  equation  (4.11) 


can  be  written  in  a  discrete  form  as 
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AW  =  favg  (Z,t)  .  At 


(4.12) 


Now,  if  we  consider  that  at  time  t  the  wave  W  is  at  the 


Kt^1  node,  then 


AW  =  WK(t)  -  WK-1  (t  -  T) 


(4.13) 


In  section  (3.6)  it  was  shown  that  the  second  formulation 

is  superior  to  the  first  formulation,  so  we  adopt  the  former  in 

approximating  f(Z,t)  I  so  that  we  take 

avg 


or 


f^Z’t)  'avg  2  ^ f K-l  ^  T‘}  + 


(4.14) 


Thus  we  have  for  the  forward  propagating  argument 


VK<t')  +  CK  1K^  VK-l('t  T)  PK-1  1K-1  ^ 


-  f  [  4  VKU)  VkS  lR(t>  + 


&  v.  <■-»*  ‘K-1 


K-l 


i .  "K-i 


iT,  ,  (t  -  T)] 


1  K-l 


(4.15) 


vK(t)  [i  +  f  +  ^k  +  T^c 

k  vlkck 


1  JK(t>  * 


VK-l(t  T)  U  -  2  C  . 1  +  [pK-l  2 

K-l 


At  RK-1 


ihK- 


l  ck-i 


'^K-l^  T) 


K=2, . . . ,N  (4.16) 


Similarly  for  the  backward  propagating  wave 
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Vv  ,  (t)  [1  +  -^ 


K-l 


2  C 


K-l 


r  ,  At  \-l 

-  [P^  +  2 - 


K-l 


VK(t  -  T)  [1  -  f5-^] 

K 


At  RK 

tpK  '  TWV'  ’  T> 


VL~C 


K  K 


K=2 , . . .N  (4.17) 


Let 


R 


*K 


=  At 


K 


/  L..C 


K  K 


(4.18) 


bK  =  At  C 


(4.19) 


Then  (4.16)  and  (4.17)  transform  to 


VK(t)  [1  +^f  ]  +  (PK  +  ]  iK  - 


Vl(t  '  T)  ^  -  2 


K-l 


]  +  [pk-1  "  “P"  1  1K-l(t  "  T)  (4'20) 


VK-l(t^  ^  +  2  ^  ^PK-1  +  2  ^  XK-1 ^ 


bK  aK 

VK  (t  -  T)  [1 - 2  1  -  [pK - 2  ]  iK(t  ”  T) 


(4.21) 


Rewrite  (4.21)  between  the  Ktb  and(K  +  t)SC  nodes. 


st 


we  have 


vK(t)  [i  +  ] 


lpK  +  -T  1  1K(t) 
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VK+l(t  -  T)  U 


bK+l  _  aK-H 

2  J  PK+1  2 


iK+1(t  -  T)  (4.22) 


further  adding  (4.20)  to  (4.22)  and  subtracting  (4.22)  from  (4.20) 
we  get 


VK(t)  =  [(1  -  )  VK_x(t  -  T)  + 

bK+l  av-_i 

(!  2  )  VK+l(t  T)  +  ^PK-1  2  ;  ^-1^  T') 

^K+l 

(fK+i  -  )  iK+J(t  -  T)]  /  [2  +  bR]  (4.23) 

bK-l 

iK(t)  =  [(1  -  ~Y~  )  vK_1(  t  -  T)  - 

av,_  I 

(i  2  ^  VK+l('t  "  +  ^PK-1  2  ^  1K-1  ^  “  T') 


+  (P 


K+l 


aK+l. 

2  ' 


XK+l(t  T)]/[2pK+aR] 


(4.24) 


CHAPTER  5 


DESCRIPTION  OF  COMPUTER  PROGRAM 

The  algorithms  discussed  in  chapters  3  and  4  were  in¬ 
corporated  into  a  Fortran  IV  computer  program.  The  operation  and 
flow  of  the  main  program  are  described  in  this  chapter. 

A  sample  computer  output  listing  is  given  to  illustrate 
the  operation  of  the  program. 

5.1  Computer  program  flow  chart 

The  flow  of  the  program  is  shown  in  fig.  (5-1)  and  the 
program  listing  is  given  in  Appendix  C. 

The  program  reads  in  data  cards,  describing  the  system 
parameters,  accuracy  limits  and  type  of  transmission  line  whose 
transient  response  is  to  be  analyzed.  The  program  then  writes  out 
all  input  parameters,  type  of  transmission  line  (uniform  or  non- 
uniform)  and  maximum  permissible  relative  error  (in  case  of  uniform 
lines ) . 

In  the  case  when  a  uniform  line  is  analyzed,  the  program 

determines  the  appropriate  number  of  sections  n.  If  the  computed 

value  of  n  is  less  than  10,  the  program  sets  n  to  be  equal  to  10 

and  the  main  loop  is  entered.  Also  if  the  computed  value  of  n  is 

less  than  1000  the  main  loop  is  entered.  If  the  computed  value  of 

n  is  larger  than  1000,  the  program  writes  out  a  message  indicating 

that  the  number  of  sections  has  exceeded  the  maximum  permissible 

value  and  the  program  stops.  It  should  be  noted  that  this  limit 
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can  be  increased  if  desired  by  a  minor  change  in  the  program  at  the 

expense  of  higher  storage  requirements. 

If  the  line  is  nonuniform,  the  calculation  of  n  is  skipped 

and  the  main  loop  is  entered.  Here  n  is  given  in  the  data. 

The  first  step  in  the  main  loop  is  to  divide  the  line’s 

physical  length  into  one  thousand  sections  of  equal  length.  It  is 

assumed : that  the  line  is  of  unit  length.  The  vector  Z  of  physical 

K 

length  between  the  sending  end  and  the  node  is  created.  Then 

the  corresponding  electrical  length  vector  Y  is  computed  using  a 

trapezoidal  rule  formulation. 

The  line’s  electrical  length  £  is  then  computed  and  the 

/  \ 

vector  of  normalized  electrical  length  is  created  Y^  by  dividing 

K 

Y  through  z, . 

K 


The  program  then  computes  the  appropriate  time  step  At 
as  the  out- come  of  dividing  the  total  electrical  length  by  the 
number  of  sections  n.  The  line  is  then  divided  into  n  sections  of 


equal  electrical  length  (delay)  by  creating  the  vector  Y  accord¬ 
ing  to  the  rule  Y^  =  At.j. 

Further  the  physical  length  vector  corresponding  to 

Y  is  computed  using  a  linear  interpolation  formula.  With  the  vector 

Z.  created,  the  vectors  p.,  Cof ,  . ,  Cof_.,  Cof_.  and  Cof ,  .  are 
J  3  lj  2j  5j  6j 

computed  as  a  prior  step  to  applying  the  current  and  voltage  relations, 
The  initial  values  of  current  and  voltage  at  each  node 


are  then  stored  as  zeros  in  the  C  and  V  arrays. 


- 


START 


\  READ  7 

\  SYSTEM  / 

V / 


NO 


COMPUTE  a  ,  n 


NO 


CREATE 
VECTORS  ZK, 
SK  AND  Yk 

r  =  Y 
T1001 

-Y. 

YES 


YES 


o 

ii 

c 

40 


CREATE  VECTOR 
NORMALIZED  YK 


COMPUTE  VOLTAGE 
AND  CURRENT  AT 
EACH  NODE  AT  TIME  t 


BOUNDARY 

CONDITIONS 


CL 


o 

O 

_ i 

uu 


5 


INDEX  LOOP 


OLD  -  NEW 
VALUES  EXCHANGE 


STORE  CRT 
VALUE  FOR 

PLOTTING 


PLOT  LD  CRT.,  INPUT 
CRT. ,  p  vs  Z  AND 


Y„  VS  Z 


FIG.  5.1  COMPUTER  PROGRAM  FLOW  CHART 


42 


In  this  program  C  and  V  denote  old  values  of  current 
and  voltage,  while  CC  and  VV  denote  new  values  of  these  variables. 

By  old  and  new  values,  it  is  meant  values  at  time  instants  (t-n.At) 
and  at  t  respectively.  The  next  step  is  to  compute  the  values  of 
new  currents  and  voltages  at  each  node  using  old  values.  Then  the 
program  replaces  the  old  values  by  new  values  to  determine  the 
currents  and  voltages  at  the  next  time  instant  and  so  on  until  the 
required  time  interval  is  finished.  Thus  the  program  destroys  old 
values  of  voltages  at  every  node  and  currents  at  every  node  except 
at  the  receiving  end. 

The  program  then  writes  out  the  values  of  currents  at  the 
sending  end,  one  quarter  line  length,  half  line  length,  three  quarter 
line  length  and  receiving  end  at  every  fifth  time  instant. 

The  last  step  in  the  program  is  to  plot  the  input  current 
wave,  load  current  wave  versus  time.  Also  the  program  will  plot 
the  variation  of  characteristic  impedance  p  and  normalized  electrical 
length  with  the  physical  length. 

5-2  User  supplied  data 

The  computer  program  user  must  supply  data  to  the  program 
describing  the  characteristics  of  the  system,  the  calculation 
accuracy  or  number  of  sections. 

These  data  are  supplied  to  the  program  by  means  of  data 


cards  at  the  end  of  the  Fortran  Card  deck. 


. 
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The  parameters  and  the  formats  on  each  data  card  are 
described  in  this  section. 

The  first  data  card  contains  the  load  resistance  RR  and 
the  source  shunt  conductance  GG.  The  format  is  2E16.8. 

The  second  data  card  contains  XNN  and  XMEW.  These  are 
the  number  of  sections  n  and  the  relative  error  allowed  y.  If  the 
value  of  y  supplied  is  zero,  the  program  takes  the  line  as  a  non- 
uniform  one.  However,  in  this  case  the  value  of  n  must  be  entered. 

The  format  is  2E16.8. 

The  third  data  card  contains  RESI,  R0VAL  and  GVAL ,  these 
values  are  to  be  entered  if  the  line  is  uniform.  RESI  is  the  value 
of  the  resistance  per  unit  length,  R0VAL  is  the  value  of  the 
characteristic  impedence  and  GVAL  is  the  value  of  the  shunt  con¬ 
ductance  of  the  line  per  unit  length.  The  format  is  3E16.8. 

The  fourth  data  card  contains  XLO,  CO,  RO ,  GO.  These 
are  the  amplitudes  of  inductance,  capacitance,  resistance  and  conduct¬ 
ance  per  unit  length  functions  describing  the  variation  of  these 
parameters  with  the  physical  length.  The  format  is  4E16.8. 

The  fifth  data  card  contains  XLE1,  XLE2 ,  XLE3 ,  AMP  and 
SLOPE.  These  parameters  describe  the  source  current  characteristic 
as  a  trapezoidal  wave.  This  wave  is  shown  in  fig.  5-2  with  the 
meaning  of  each  of  the  above  mentioned  parameters. 

It  is  noted  that  SLOPE  =  AMP/XLE1.  The  format  is  5E12.-4. 
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FIG.  5.2  TRAPEZOIDAL  WAVE  FORM 

In  addition  to  the  above  mentioned  five  data  cards  the 
following  functions  should  be  supplied. 

The  first  function  is  FL(X),  which  describes  the  variation 
of  the  line’s  inductance  with  the  length  X.  The  second  function 
is  FC(X)  for  the  capacitance.  The  third  function  is  FR(X)  for  the 
resistance.  The  fourth  function  is  FG(X)  for  the  shunt  conductance, 
and  the  fifth  function  is  the  source  current  function  SOUR(X). 

As  an  example  the  following  functions  were  assumed  for  the 
above  mentioned  variations. 

FL (X)  =  XLO*  (2.0  +  SIN(ttX)  ) 

FC (X)  =  CO  *  EXP  (-X) 

FR(X)  =  RO  *  Log  (1  +  X) 


FG  (X)  =  GO 
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5-3  Numerical  study  of  the  effect  of  n  on  results 

The  program  described  in  this  chapter  was  equipped  with 

a  facility  to  compute  and  plot  the  transient  response  of  the  given 

line  for  three  different  values  of  n. 

This  is  facilitated  by  the  introduction  of  a  three  step 

D0  loop.  The  first  for  the  given  value  of  n  say  nQ,  the  second 

for  n  =  5n  and  the  third  for  n  =  20n  . 

o  o 

Another  trial  was  for  n  ,  1.5n  and  2n  respectively. 

o  o  o 

The  computed  values  of  currents  at  the  locations  on  the 
line  previously  mentioned  and  plots  of  load  current  versus  time 
are  shown  in  figs.  (5-3)  and  (5-4). 
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CONCLUSION 

A  method  for  the  evaluation  of  the  transient  response  of  a 
nonuniform  lossy  transmission  line  has  been  presented  in  this  thesis. 
The  pair  of  partial  differential  equations  which  describes  the 
voltage  and  current  relations  on  the  line  is  transformed  into  a  pair 
of  ordinary  differential  equations  on  two  characteristic  curves 
using  the  method  of  characteristics. 

A  stepped  line  approximation  is  used  to  analyze  the 
transient  response  of  the  given  nonuniform  line.  The  concept  of 
electrical  length  is  employed  in  dividing  the  line  into  a  number  of 
equal  delay  sections.  The  set  of  difference  equations  describing 
the  set  stepped  line  is  suitable  for  digital  computer  solution. 

The  main  advantage  of  this  method  is  that  numerical 
techniques  such  as  the  Runge-Kutta  method  are  entirely  avoided. 

The  computational  procedure  involving  the  use  of  a  digital 
computer  is  illustrated  for  a  specific  distributions  of  L(x),  C(x), 
r (x)  and  g(x). 

There  is  still  much  work  which  has  to  be  done  in  this  area. 

For  instance,  it  will  be  worth  while  to  investigate  the 
possibility  of  obtaining  an  estimate  of  the  proper  number  of  sections 
to  achieve  a  given  accuracy  in  the  case  of  nonuniform  lines. 
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APPENDIX  A 

BASIC  EQUATIONS  FOR  A  UNIFORM  LOSSLESS  LINE 

Consider  a  uniform  lossless  transmission  line  of  length  1 
as  shown  in  figure  (A.l).  The  following  two  equations  hold  true  for 
the  voltage  and  current  values  w.r.t.  the  position  on  the  line  x 
and  time. 


3  V(x,s) 
3  x2 


=  s  L  C  V  (x,s) 


(A-l) 


9  -(-y^-S)  =  S2  L  C  I  (x,S> 

3  x 


(A-2) 


Where  L  and  C  are  the  inductance  and  capacitance  per  unit  length  of 
the  line  respectively,  and  s  being  the  Laplace  operator. 

The  solution  of  A-l  and  A-2  takes  the  form 


-sx 
v 


V(x,s)  =  F(s)  e  +  G(s)  e 


sx 

v 


(A-3) 


where 


-sx 
v 


p  I(x,s)  =  F(s)  e  -  G(s)  e 


sx 

v 


v  = 


(LC)  2 

is  the  propagation  velocity. 
P  =  (L/C)^ 


(A-4) 


(A-5) 


(A-6) 


is  the  characteristic  impedance  of  line. 

Applying  A-3  and  A-4  at  the  boundaries,  x=0  and  x=l  we  obtain 

-SL  -sL 

V(l,s)  +  pi  (l,s)  =  V(0,s)  eV  +  pi  (0  ,s)  eV 


(A-7) 


-sL 
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V(0,s) 


Pi  (0,s) 


-sL 

V(1 , s)  eV 


pl(l,s)  eV 


Letting 


v 


Then  we  get  in  the  time  domain 


(A-8) 


V(l,t)  +  pi(l,t)  =  V(0,  t-T)  +  pi(0,  t-T)  (A-9) 

V(0,t)  -  pi (0 , t)  =  V(l,  t-T)  -  pi(l,  t-T)  (A-10) 

It  should  be  noted  that  the  classical  transmission  line 
analysis  adopted  here  assumes  the  following: 

1.  The  transverse  electric  field  components  in  the  conductor 

are  negligible  compared  with  the  axial.  This  is  equivalent  to  the  assumption 
that  displacement  currents  in  the  conductor  are  negligible  compared  to 

conduction  currents.  This  is  stated  mathematically  as  —  >>1. 

toe 

2.  The  axial  electric  field  components  in  the  dielectric  are 

small  compared  with  the  transverse.  This  is  equivalent  to  the  assumption 

that  the  characteristic  impedance  of  the  dielectric  is  much  greater  than 

R 

g 

the  skin  effect  surface  resistivity  of  the  conductor  or  —  <<1. 

n 

These  inequalities  are  nearly  always  satisfied  by  the  materials 
of  common  transmission  lines,  but  if  they  are  not,  one  must  examine 
critically  any  results  predicted  by  the  usual  transmission  line  equations. 
Also,  the  longer  the  duration  of  the  input  current  function,  the  better 
these  inequalities  are  met. 
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-1  A  Uniform  Losseless  Transmission  Line 
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APPENDIX  B 

DERIVATION  OF  EQUATION  (3-50) 

Consider  an  incremental  length  Ax  at  a  distance  x  from 
the  sending  end  of  a  uniform  lossy  transmission  line  with  length  h 
This  incremental  length  is  shown  in  fig.  B-l,  then  we  can  write 
the  following  equations  : 


=  (sL  +  R)  I 

3x 


(B-l) 


-fi  -  (SC  +G)  V 


(B-2) 


where  R,  L,  C  and  G  are  resistance,  inductance,  capacitance 
and  conductance  per  unit  length  respectively,  s  is  the  Laplace 


operator , 


yields  : 


Partial  differentiation  of  (B-l)  and  (B-2)  w.r.t.  x 


32V 

3x2 


=  (sL  +  R) 


31 

3x 


(B-3) 


3 2l  .  ,  -  .  rs  3 V 

~2  ‘  +  G)  3x 

3x 


(B-4) 


c\  t 

Substituting  for  —  and  —  in  (B-3)  and  (B-4)  from  (B-2) 

dX  o  X 


and  (B-l) 
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respectively  yields 


2 

3  V 


9x 


2  =  (sL  +  R)  (sC  +  G)  V 


(B-5) 


92I 


9x 


2  =  (sL  +  R)  (sC  +  G)  I 


(B-6  ) 


Let 


Y  (s)  =  (sL  +  R)  (SC  +  G) 

Then  (B-5)  and  (B-6)  can  be  rewritten  as 


(B— 7 ) 


92V  _  2  .  . 

- 2  =  Y  (s)  V 

9x 


(B-8) 


2 

3  1  2  f  .  _ 

- 2  =  Y  (s)  I 

9x 


(B-9 ) 


The  solution  to  (B-8)  is 


V(x,s)  =  F(s)  e  YX  +  H(s)  eYX 


(B-10 ) 


where  F(s)  and  H(s)  are  functions  to  be  determined  such  that 
boundary  conditions  along  the  line  are  satisfied.  Using  (B-10) 
and  (B-l)  we  get  for  I(x,s)  : 


!(x,s)  =  [F (s)  e  YX  -  H(s)  eYX] 


(B-ll ) 


In  (B-10)  and  (B-ll),  forward  propagating  components  of 


voltage  and  current  are  the  ones  having  negative  exponents,  and 


, 
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backward  propagating  components  are  those  with  positive  exponents. 

Thus  if  suffix  f  denotes  forward  components  then  we  have 

Vf(x,s)  =  F (s)  e  YX  (B-12) 

VX’S)  =  sL+R  F(s)  e’YX  (B-13) 

now  at  x  =  0  we  have 


Vf (0 , s)  =  F(s) 

V°>s)  F<s> 

Thus  (B-12)  and  (B-13)  rewritten  as 
Vf(x,s)  =  Vf(0,s)  e  YX 

If(x,s)  =  If  (0 ,  s)  e  YX 

Now  by  definition  of  W(x,s)  we  have 

W(x,s)  =  V  (x,s)  +  pi  (x,s) 
f  t 

then  using  (B-16)  and  (B-17)  one  gets 

W(x,s)  =  W (0  , s )  e“YX 

Let 

WR  =  W  (0  ,  s  ) 


(B-14) 

(B-15) 


(B-16) 


(B-17) 


( B-18) 


(B-19)' 


(  B-20) 


\+l=  W(h»s) 


(  B-21) 


- 

I 
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then 

WK+1  =  WK  6XP  ["Yh]  (B-22) 

Now  we  have  by  (B-7) 


Y  (s)  = 

[s2LC  +  (RC  +  GL)  s  +  RG]^  (B-23  ) 

then 


Y  (s)  = 

s/Zc  [1  +  i  (  |  +  |)  +  5|_ ]**  CB-24) 

s  5  LC 

If 

the  time  t  is  small,  then  s  is  very  large  ( s>°°  as 

t  ->  0) ,  then  we  can  write 

Y(s)  =  sfc  [1  +  i  (  5.  +  I)]®5  (B-23 

Using  the  binomial  theorem  one  gets 


Y  (s)  = 

sjLC  [1  +  yr  (  f-  +  |)]  (  B-26  ) 

or 

y(s)  =3^  +  |  [  J+G  p]  (  B-27  ) 


h. Y  (s) 

=  s.h.  J~L C  +  4  [  —  +  G  h.p  ]  (  B-28  ) 

2  p 

Now  by  definition  of  At  as  the  electrical  length  of  the  section 
with  physical  length  h  and  letting 

r  =  R.h  (  B-29  ) 


(  B-30  ^ 


g  =  G.h 


we  have 


h.y(s)  =  At.  s  +  +  g  p] 


(B-3l) 


Thus 


WK+1  =  WK  exp  [  -  At-S  "  \  (  f  +  S  P)1 


(B-32 ) 


which  is  equation  (3-50) 
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FIG.  B-l  INCREMENTAL  LENGTH  OF  A  TRANSMISSION  LINE 


APPENDIX  C 


PROGRAM  LISTING 


I 


0002 

0003 

0004 

0005 

0  005 

0  C  0  7 

0  OOP 

0  C0<? 

0  0  1  c 

0  0  11 

0  0  12 

0  0  13 

0  0  14 

0  0  15 

o  o  i  e 

0  0  17 

0018 

0  0  1  s 

0020 

0  02  1 

0022 

0  023 

0  024 

0025 

0026 

002P 

0029 

0030 

003  1 

0  032 

0033 

0034 

0035 

0  C  3  7 

0  0  39 

0  04  0 

0  04  1 

0042 

004  3 

0  04  4 

0045 

0  04  6 

0  C4  7 
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c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


1 ....  1 ....  1 

GENERAL  TRANSMISSION  LINE  *4******************* 
TRANSIENT  ANALYSIS  ************************* 


C I K  )  . VI K )  . . .OLD  CURRENTS  AND  VOLTAGES  AT  THE  SECTION  BOUNDARIES 
CC( K )  , VVI K  )  . . .NEW  VALUES  CALCULATED  FROM  C < K ) . V ( K  > 

LINE  DIVIDED  INTO  NK-l  SECTIONS 

TIMENSICN  Y  AB  (  10 04  ). XI  1<  1004), YBB I  1004), X22I  1004), YCB1  1004)  ,  X33(  10 


204  > 

Cl  MENS  I  ON  C I  1 000 ) ,CC(  1000),V(500),VVI500),X(I002,4).WORK( 
DIMENSION  F  Y (  1001  )  ,  Z(  1001  >.DEL<  100 l  ).SY(  1001  )  ,  S7(  1002.3) 
DIMENSION  PU( 1001 ) ,YN( 1002) ,Z1 ( 1002) 

DIMENSION  74(1002), COE1J501  ),COF2(501  ), COF5( 501  ).C OF  6(501 
DIMENSION  P  CO  (  1002.4  )  ,Y (  1 002,4)  «Y3(  1002.4)  ,Y4(  1002  ) 
DIMENSION  RAB(I04).RBB(504),RCB(1004),Z11(104), 722(504). Z 
EXTERNAL  FL , E C . FR . F G , SOUR 
COMMON/BCB/XLO.CO.RO, GO . XLE 


1  024  ) 


) 

3  3 (  1  004  ) 


COMMON /SOS/ XLE 1 . XLE 2 , X LE 3 , A MP , SL OPE 
WORK ( N ) .  .  .  . AUX IL I  ARY  VECTOR  FOR  PLOTTING 
CALL  PLOTS  (WORM  1  )  ,4096) 

CALL  PLOT  (5. 0.5. 0,-3) 

OUTPUT  RESISTIVE  LOAO=FR 
PE  AD (  5.901  )  RR.GG 
901  FORMAT ( 2F 1 6.8 ) 


PEAO( 5,90? ) XNN, XMEW 


902  FORMAT ( 2E l 6 . 8 ) 

READ(  5.903)  R  F  S  I  ,  RO  V  AL  .  C,  V  AL 

903  FORMAT! 3E 1 6 .8 ) 

READ! 5.904)  XLO.CO.RO.GO 

904  FORMAT ( 4E 16. P ) 

WRITE(6,910)  XLO.CO.RO.GO 

910  FORMAT! 4E 10.4) 

READ! 5.905)  XLE 1 . XLE 2 , XLE 3 , AMP . SLOPE 

905  FORMAT! 5E 12.4 ) 

IF  (XMEW. EQ. 0.0)  GO  TO  666 

ALF A=(  ( RE SI /ROVAL  )  «• ( GV AL *ROV AL  )  )/2. 

XNN  =  ALFA* SORT (  ( ALFA  )  / (  1 2 . *XMFW)  ) 
WRITF(6,906)  RFSI.ROVAL.GVAL. XMEW. ALF A , XNN 

906  FORMAT (  1 X  ,6E 1 6 . 8 ) 


NN- XNN 

WPITE(6,907)  NN 
907  FORMAT ( 30  X . 1  IP) 

IF  (NN.LT.10)  NN=10 
IF  ( NN . GT  .10  0  0)  GO  TC  777 
666  F  Y (  1  )  =  0.00 

CO  1  K  1  =2 . 1  CO  1 

DEL (  1  )  —  SCOT ((FL(O.O)  )*(FC(0.0>)> 

XK=K  1 

7 (K 1 >-( XK- 1 . ) *0.00 1 
Z  K  =  Z  (  K  1  ) 

DFL ( K  1  )=SORT((FL(ZK)  >*(FC(7K)  >  ) 

1  FY(K1 )=FYIK1-1 ) *0 .5* ( DEL ( K 1 ) FDEL (K 1- 1 ) ) 
T  OU-F Y (  10  01  ) -F  Y (  1  ) 


. 


0048 

0  C  4  5 

0050 

005  1 

0052 

0053 

0  054 

0055 

0057 

0  0  5  <3 

0060 

0  06  t 

0  062 

0063 

0  C  6  4 

0065 

0066 

0067 

0  0  6  9 

0  07  1 

0073 

0075 

0  C  7  6 

0077 

0076 

0  0  7  9 

0C8C 

008  1 

0082 

0  0  83 

0084 

0065 

O  086 

ooe7 

0088 

0089 

0090 

009  1 

0  092 

0093 

0  094 

0  C  9  5 

0  096 

0097 

0  09  8 

0  C  9  9 

0  1  0  C 

0  10  1 

0  102 

0  103 

0  104 

0  106 

0  106 
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DO  2  I X  =  1  ,  1 000 

1X1=IX>1 

2  YN(  I  X  )  =  FY (  I  X 1  l/TOU 
DO  3  I Z=1  .  1 000 

17 1= I Z+l 

3  71  1  I  7 )  =  7 (  171  ) 

DO  7  INDEX=1,3 

IF  (  INDEX  .EQ.  2  )  XNN=l.S*XNN 
IF  (  INDEX  .EO.  3  )  XNN=(4.0*XNN>/3.0 
10  DELT=TOU/XNN 
NN=XNN 
N P=NN— 1 
DO  4  J 1  =  1  ,NP 
5  Z  <  1  »  INDEX  >  =  0.0 
X  J  =  J  1 

SY(  J 1  )  =  X J  *DELT 
CO  4  K2= 1 ,1000 


IF1SY1J1  ) 

.  GT  . 

FYIK2+1 ) )  GO  TO 

4 

I F ( SY ( J 1 > 

•  LT  . 

FYIK2+1)  .AND. 

SY ( J1 > .GT . 

FY(K2) 

) 

Y1 =FY( K2 ) 

IFISYIJl  > 

.LT  . 

FYIK2+1  )  .AND. 

S  Y ( J 1  >.GT. 

F  Y  (  K  2  ) 

> 

Y6=F  Y( K.2  + 1  > 

IF(SY(Jl) 

•  LT  • 

FY1K2+1 )  .AND. 

S Y ( J 1  ) .GT  . 

F  Y  (  K  2  ) 

> 

KNN-K2 

J1  1  1  =  J1  +  1 

SZUUl  ,INDEX)  =  Z(KNN)  +  (0.001*((SY(J1)-Y1  )/(Y6-Yl)  >) 

4  CONTINUF 

SZ (MN4  1 . I NDEX )= 1 .0 

FX=NN+1 

NY=NN*2 

CO  5  J5-1 .NX 

XZ J  =  SZ (  J5  .  I NDEX  > 

XL  =  FL ( X 7  J  ) 

XC=FC(XZJ  > 

XR=FR(XZJ  > 

XG=FG( XZ J  > 

PO( J5  >  =  SQRT ( XL/XC  > 

COF  l< J5>=1.-(DELT4XG)/(2.*XC> 

C  OF  2 ( J5 )=RP( J5  >-( DELT4XR >/<  2.  *( SORT! XL  4  XC )  )  ) 

C0F5( J5  >= (2  .  >  4< 1  . ♦ (DELT4XG  )/(2.*XC >  ) 

COF6(J5)=(2.)4(RO(J5)+(DELT4XR)/(2 . 4( SORT ( XL4XC ) ) ) > 

5  FOO ( J5. INOE X >=RO< J5 > 

c  444  INITIAL  CCNDITICNS  CN  TRANSMISSION  LINE  444 
DO  6  K  3= 1  ,NN 
C (K3)  =  0 .0 

6  V(K3)=0.0 

c  444  ITERATIVE  CALCULATION,  L=DISCRETE  TIMF  *44 
NH=5*NN 
CC  7  L= 1 . NH 
XI E=L/XNN 
4  p  =  NN- 1 
DO  8  J=2.NP 

A=  V  (  J—  1  >  4 COF  1  (  J-  1  )  +  C  (  J-  1  ).*COF2<  J-l  > 

P=V( J4l  )*COFl  ( J4  1  >-C< J+l  )*COF2< J+l  ) 

V V( J ) =( A+P ) /C0F5 ( J ) 

CC  <  J  >  =  (  A-P  > / C 0F6 ( J  > 

P  CONTINUE 

C  44*  POUNDARY  CONDITIONS  *44 

XA=0.5*(C0F5(  I  )  +GG4CCF 6 (  l  )  ) 


I 
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I  SN 

0  1  07 

VV(1)  =  (0.54{  SOUR  (  XLE)AC0F6<1>)4(V<2)AC0F1<2))-(C(2>*C0F2(2))>/XA 

I  SN 

o  i  c  e 

CC< 1 l-SOUP(*LE)-GG»VV( 1 1 

I  SN 

0  10° 

XB=0.54(COF5(NN)4(COF6(NN)/PR) ) 

I  SN 

0  110 

VV(NN)=( ( V ( NP ) 4COF 1 ( NP ) ) 4 ( C ( NP ) 4COF21 NP ) ) )/XB 

l  SN 

0  111 

CC  (  NN)  =  VV  (  KM  /RR 

c 

*44 

OLO  -  NEW  VALUES  EXCHANGE  444 

I  SN 

Oil? 

CO  9  M= 1 , NN 

I  SN 

0  113 

C I M) =CC ( M ) 

I  SN 

0  114 

9 

V( M l^VVIM  ) 

I  SN 

0  11S 

444 

X ( L •  I NOEX  )  =  XLE 

I  SN 

out 

Y  <  L  .  I NDEX )=C ( NN) 

I  SN 

0  117 

Y3(L ,  INDt  X  )  =  SOUR (  XLE  ) 

c 

PRINT  THF  CURRENTS  AT  THE  EVERY  FIFTH  TIME  INSTANT: 

1  SN) 

0  1  1  P 

IF(L/545.E0.L)  WRITE  (6.908)  C (  1  )  , C ( NN/ 4 )  . C ( NN/ 2 ) . C ( 3* NN/ 4  )  ,  C ( NN ) 

1  SN 

0  120 

008 

FORMAT  (  1  X  .5F20.7  ) 

1  SN 

0  1  ?  1 

7 

CONTINUE 

I  SN 

0  1  22 

RO I 1=0 . 0 

I  SN 

0  123 

YBI 1=0.0 

I  SN 

0  124 

N AX  =  ( 54NN  )/? 

I  SN 

0  1  25 

CO  71  L= 1 . NAX 

I  SN 

012t 

YAB(L)=Y(L . 1 ) 

I  SN 

0  127 

YBI  1  =  A  P  A  X 1  <  Y  A  E ( L )  .YBI1  ) 

I  SN 

o  i  2e 

71 

X  1  1  ( L  )  =  X ( L ,  1  ) 

1  SN 

0  120 

1  CN  =  NN/2 

I  SN 

0  130 

00  75  L-l . I CN 

I  SN 

0  13  1 

P AB(L )  =  ROC (Li  1  ) 

I  SN 

0  1  32 

RBI  l  =  AMAXl  (RAE'(L).RBIl  ) 

I  SN 

0  133 

7 1  l(L)=SZ(L  .1  ) 

I  SN 

0  134 
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C  C  NT INUE 

I  SN 

0  135 

YCI 1=0.0 

I  SN 

0  135 

PC  I  1  =  0 . 0 

I  SN 

0  137 

NBX= ( 1 54NN )/4 

T  SN 

0  1  3  P 

DO  72  L= 1  .NBX 

ISN 

0  1  30 

YPO(L  )=Y(L  .2) 

I  SN 

0  140 

YC I  1 =AMAX 1  (YPB(L).YCIl  ) 

l  SN 

0  14  1 

72 

X22(L)=X(L .2) 

1  SN 

0  14? 

JCN=( 34NN 1/4 

I  SN 

0  143 

CC  77  L=t  ,JCN 

l  SN 

0  14  4 

PPfl(L)  =R00  (  L  .2  ) 

I  SN 

0  145 

PC  11  =  AMA  X 1  (RBB(L).RCIl) 

1  SN 

0  145 

722IL )=SZ (L  ,2  ) 

ISN 

0  147 

77 

CUNT INUE 

I  SN 

0  1  4  P 

Y  D I  1  =  0.0 

I  SN 

0  140 

PO  I  1 =0 .0 

1  SN 

0  150 

N  C  X  =  5*NN 

I  SN 

0  15  1 

CO  73  L  = 1  . N  C  X 

I  SN 

0  152 

YC  B ( L ) = Y ( L .3) 

I  SN 

0  15  3 

YDI1=AMAX1(YCB(L).Y0I1) 

I  SN 

0  154 

73 

X  3  3  (  L  )  =  X  (  L  .  3  ) 

I  SN 

0  155 

K  0N  =  NN 

I  SN 

0  1  55 

CO  78  L  = 1  ,KCN 

1  SN 

0  157 

PC  B ( L )  =  PCO( L . 3 ) 

I  SN 

0  1  5P 

EDI  1  =  A  M  A  X 1  (RCB(L)  .RDI  1  ) 

I  SN 

0  1  50 

? 13(1 )=SZ(L.3> 

1  SN 

0  15C 

76 

CONT I NUE 

I  SN 

0  16  1 

yo  !G= AMAX 1  ( YB 1 1  .YCI  1  . YD  I  1  ) 

- 


0  162 

0  163 

0  164 

0  165 

0  166 

0  167 

0  168 

o  i  es 

0  170 

0  17  1 

0  172 

0  l  73 

0  1  74 

0  175 

0  176 

0  177 

0178 

0  1  79 

0  180 

o  i  e  i 

0  1  82 

0  1  83 

o  i  e4 

o  i  e? 

0  1  86 

0167 

0  188 

o  i  es 

0  1  90 

0  19  1 

0  192 

0  1  93 

0  1  94 

0  195 

0  196 

0  197 

0  198 

0  1  99 

0200 

020  1 

0  202 

0203 

0204 

0  205 

0206 

0207 

o  2oe 

0209 

02  1  C 

02  11 

02  1  2 

02  13 

02  14 

0  2  l  K 

02  16 

0217 

02  1  8 

0219 

0220 

0  22  1 

0  222 

0223 

0224 
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RBIG=AMAX  1  C RB  I  1  .RCI  1  ,RDl 1  > 

WR1TFC6.74)  YBIGt RB I G 
74  FORMAT(30X,2E16.8./) 

C  ♦**  PLOTTING  INSTRUCTIONS  *** 

RAR( I0N4)  1  =  0 .0 
FAB<  I  ON +  2 )=RBIG/4. 

7  1  1  (  I ON  4 1  )  =  0.0 
7  111  I0N42  1=0.2 

C  ALL  AXtSC0.0,0.0,'RC',2,4.0,90.0.RAB(ION4l),RARCION42),20.0) 

CALL  AXISIO. 0.0.0, 'ACTUAL  LE NTH*. -13.5. O.O.O.ZIICION411.Z11C ION4  2 ) 
1.20.0) 

CALL  L I NE ( Z 1  1  ,RAB  .  ICN, 1  .0 .3 ) 

CALL  PLOT (0.0. 0.0, -3 1 
PBB<  J0N41  1  =  0 .0 
RPB<  JON 4  2  )  =  RB I G/4  . 

7  221 J0N4 1  1  =  0.0 
Z  2  2 ( JON42 1  =  0.2 

CALL  LINECZ22.RBB, JON, 1 . 0, 3) 

CALL  PLOT  CO  .0 , 0. 0.-3 1 
RCBCK0N4I  1  =  0.0 
PCB(KON42)=RBIG/4. 

7  33 ( K0N4 1  1=0.0 
733<KON42)=0.2 

CALL  LlNE(Z33,RCB,KON. 1,0,3) 

CALL  PLOT (0 .0  .5.5 ,-3 ) 

CALL  SCALE  { YN  .  4 . 0 , 1 000 , 1 . 20 . 0 1 

CALL  AXIS  ( 0 .0 ,0 .0 .' YN(  I  X)  •  ,6, 4.0.90. 0 . YN(  1001  1  , YNC  1  002 1 . 20.0 1 
CALL  SCALE (  Z1  .5.0 . 1 000,  1 . 20 .0  1 

CALL  AXIS  ( 0.0 .0. 0.  'ACTUAL  L E NGTH •  , - 1 3 . 5 . 0 . 0 . 0 . Z l  < 1 0 0  l  1  ,  Z  1  (  l  0 02  1 , 2 
10.0) 

CALL  LINE  <  Zl  .VN. 1 000. 1 .0.3) 

CALL  PLOT  ( 12. 0.0. 0.-3) 

YABCNAX41  1  =  0.0 
YABCNAX42  1  =  YB  I  G/4  . 

X  I  1  (NAX41  1  =  0.0 

X I  1 ( NAX  4?  )=1  .0 

CALL  AXIS  ( 0 . 0 .0. 0 , • LOAO  C URRENT • . 1 2 . 4 . 0 , 9 0 . 0 . Y AB ( N A X 4 1 ) , Y ARC NA X 42 
11.20.0) 

CALL  AXISCO. 0.0.0,' TI WE  •  • —4  ,5.0, 0.0, XI  1  (NAX41  )  ,  X 1  1 ( NAX42 1 , 20. 0  I 
CALL  L INF C X 1 1 , YAB .NAX, 1 . 0 , 3 1 
CALL  PLOT C 0 .0 . 0.0, -3 ) 

YBBCNBX4 1 1=0.0 
YPBCNBX42  1  =  YB  I  G/4  . 

X22C  NRX4 1  )  =  0.0 
X  2 2 ( NBX 4 2 1  =  1.0 

CALL  L  INF C X22  .  YBB .NBX,  1 .0 . 3  1 
CALL  PLOT  CO  .0 .0.0 , -3 ) 

YCBCNCX41  1  =  0.0 
YCB  CNCX42  1  =  YB I G/4  . 

X33CNCX4  1  1  =  0.0 
X33(NCX42)=1 .0 

CALL  L INF < X33  .  YCB ,NCX , 1 . 0, 3 1 
CALL  PLOT  CO. 0.5. 5, -3) 

00  600  1=1.1000 
Y  4 (  I ) =Y  3  <  1.3) 

600  CCNTINUF 

CALL  SCALE (Y4 .4.0, l 000, 1 ,20.6) 

CALL  AXISIO.  0,0.0.  •  INPUT  CURRENT  13.4.0.90.  '>,YA(  |PP1  ),Y4(  1002)  .20 
1.0) 

CALL  AXISCO. 0,0.0,'  T I  FE •  . -4  ,5 .0 .0. O  ,X TI C NCX4 I  1.X33CNCX 42  1.20.0) 
CALL  I  INt ( X33 , Y4 . 1 OOC, 1 , 0, 3) 

CALL  PLOT  C 1 0 .0  .7. 0 .-3 1 
CALL  PLOTCO. 0.0. 0.999) 

CO  TO  888 

777  »P I TF ( 6 .909 1  NN 

509  FORMAT ( 30X , 49H  FUR  THESE  VALUES  OF  SYSTEM  PARAMETERS  WE  GET  NN= , 1 8 
1/  ) 

see  stop 

ENO 


. 


, 

0  C  0  2 

0  C  0  3 

0004 

0  005 

0  C  0  6 

0  007 

0002 

0003 

0004 

0  005 

0  C06 

0007 

0  C  02 

0  C  03 

0004 

0  005 

0  C  Cf 

0  00  7 

oco? 

0  C  0  3 

0  004 

0  005 

0  COf 

0007 

0  C02 

0  00  3 

0  004 

0  005 

0  007 

0  OOP 

001  l 

0013 

0  0  14 
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FUNCTION  F  C (DAN) 
CCMMCN/BCP/XLO,C.O,RO.GO,XLF 
COMMON/ SOS /XLF 1  .XLE2.XLE3.AMP, SLOPE 
FC=CO*EXP (-CAN ) 

FETURN 

END 


FUNCTION  FL  ( CAN  1 
C CM MON/BO FV  XL 0 . C 0 , RO . GO . X LE 
C  CMMON/SPS/ XLF1.XLE2.XLE3.AMP, SLOPE 
FL=XL0*(2.P+SIN(22 .0*0 AN/ 7 . 0 1 ) 
PFTURN 
F  NO 


FUNCTION  F  F  (DAN) 
f  CMMCN/BCB/XLO , CO , RO  < GO  «  XLF 
CCMMON/SCS/ XL  E 1  ,XLF2.XLE3,AMP, SLOPF 
FR=RO*ALOG< 1 .O+DAN) 

F  F  TURN 
END 


FUNCTION  FC-(OAN) 

C  OMMON/BCH/XL  0  ,  C  0  ,  RO  ,  0-0  ,  X  L  F 
CCWON/SCS/  XLF  1  ,  XL  F2  .  XLE3  .  AMP  ,  SLOPE 
F  G=  GO 
F  F  T  U  P  N 
END 


FUNCTION  SOL'R(XAX) 

COMMON/ FI  CU/XL  0  iC  0  »R0  .GO  •  XLF 
COM MON/ SOS/ XL El . XLE2, XLE 3 . AMP , SLOPE 
!F< XAX.LT  .XLE1  >  SOUR  =  S L OPE *X A X 
IF(XAX.LT.XLE2.AND.XAX.GT.XLE1)  SOUR- amp 

IF ( X  AX . LT  .XLE 3 . AND. XAX . GT .XL  F2 )  SOUR= AMR*<  (  XL E 3- X A  X )  / (  X LF  3- XL E2 )  ) 

IF< XAX.GT  .XLE3 )  SCUR=0.0 

PFTURN 

FND 


